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Abstract 



We have adapted the anelastic spectral code of iBarranco fc Marcus! (120061 ) 
to simulate a turbulent convective layer with the intention of studying the 
effectiveness of turbulent eddies in dissipating external shear (e.g. tides). 
We derive the anelastic equations, show the time integration scheme we use 
to evolve these equations and present the tests we ran to confirm that our 
code does what we expect. Further we apply a perturbative approach to find 
an approximate scaling of the effective eddy viscosity with frequency, and 
find that it is in general agreement with an estimate obtained by applying 
the same procedure to a realistic simulation of the upper layers of the solar 
convective zone. 
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1. Introduction 

Dissipation of stellar tides and oscillations is often considered to be mainly 
due to the turbulent flow in their convective zones. Usually the effects of the 
turbulent flow are parametrized by some sort of effective viscosity coefficient. 
Clearly the situation is not as simple as that and the usual "fix" is to allow 
this viscosity coefficient to depend on the perturbation being dissipated, most 
notably its frequency and perhaps direction of the shear it creates. 

Completely analytical treatments start by assuming a Kolmogorov spec- 
trum for the turbulent flow and combine it with some prescription for the 
effectiveness of eddies in dissipating perturbations of the given period. Since 
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Kolmogorov turbulence is isotropic the direction of shear is unimportant in 
such prescriptions. Two s uch prescrip t ions h ave been used. 

The first, proposed by IZahnI (jl966l . Il989l ). states that when the period of 
the perturbation (T) is shorter than the turnover times (r) of some eddies, 
their dissipation efficiency should be proportional to the fraction of a churn 
they manage to complete in half a perturbation period, in other words, the 
effective viscosity coefficient scales like: 



mm 



2r) 
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Th e second prescription is due to lGoldreich &: Nicholson! (119771 ) and lGoldreich Sz Keeley 
(I1977I ). They argue that eddies with turnover times much bigger than the 



period of the perturbation will not contribute appreciably to the dissipa- 
tion, and hence the effective viscosity should be dominated by the largest 
eddies with turnover times r < T/27r. Then the Kolmogorov prescription of 
turbulence predicts that the effective viscosity will scale as: 
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Zahn's prescription has been test ed against tidal cir c ulariz ation times 
for binaries containing a giant star (jVerbunt fc Phinneyl . 119951 ). and is in 
general agreement with observations. Also Zahn's prescription is in bet- 
ter agreement with observ ed tidal dissipation of binary stars in clusters 
(IMeibom &: Mathied . 120051) and with the location of the red edge of the 
Cepheid instability strip (iGonczil . Il982l ). 



The less efficient prescrip t ion ha s bee n used successful l y bvlGoldreich &: Keeley 



(119771 ). iGoldreich fc Kumad (Il988f ) and iGoldreich et all (Il994f ) to develop a 
theory for the damping of the solar p-modes. In this case the more effective 
dissipation would require dramatic changes in the excitation mechanism in 
order to e xplain the observed amp htudes. 

Finally iGoodman fc OhI (119971 ) developed a perturbative derivation of the 



convective viscosity, which for a Kolmogorov scaling, gives a result that is 
closer to the less efficient Goldreich & Nicholson viscosity than it is to Zahn's. 
While providing a firmer theoretical basis for the former scaling, this does 
not resolve the observational problem of insufficient tidal dissipation. 

The development of 2D and 3D simulations of solar convection hint at 
a possible resolution of this problem. The convective flow that these simu- 
lations predict is fundamentally very different from the assumed Kolmogorov 
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turbulence (ISofia fc Chanl . ll984j : lstein fc Nordlundll989l : lMalagoli et al.l . ll990h . 



The first major difference is that the frequency power spectrum of the ve- 
locity field is much flatter than the Kolmogorov power spectrum, and hence 
one might expect that the dissipation will decrease significantly slower as fre- 
quency increases relative to the Kolmogorov case. Another major difference 
is that the velocity field is no longer isotropic and hence one would expect 
it to react differently to shear in different directions. That is, if we would 
use an effective viscosity coefficient, it should be a tensor and not a scalar 
quantity. 



As a first step in investigating that possibi htv (IPenev et al. 



20071), we 



applied the perturbative approach developed by iGoodman fc OhI ( 1l997l ) to a 
numerical model of solar convection (Robinson et al. 2003) to find the scaling 
of the components of the effective viscosity tensor with frequency. Somewhat 
unexpectedly, we found that the scaling closely follows Zahn's prescription, 
even though when the same approach is applied to Kolmogorov turbulence 
it gives results closer to those of Goldreich and collaborators. 

That left a lot of questions unanswered. For example, are the effects of 
turbulence anything at all like that of molecular viscosity, and is taking only 
the lowest order term in a power series expansion of the energy dissipation 
rate an acceptable approximation? The perturbative approach also assumes 
that the spatial scale of the external shear being dissipated is large compared 
to all convective scales, and hence can be assumed a linear function of posi- 
tion, but it would be interesting to see how quickly dissipation efficiency is 
lost as the spatial scale of the external shear decreases. 



To tackle these questions we adapted the anelastic spectral code of lBarranco &: Marcus 



( 120061 ) to simulate a convectively unstable box in which we are able to place 
external, time dependent shear as part of the evolution equations and observe 
the effects of the turbulent convective flow directly. 

In order to make the code suitable to simulate convection we added heat 
diffusion, that allows for the supply the heat that will drive the convective 
motions. This in turn necessitated that we allow for a more general back- 
ground state slightly modified treatment of the pressure and the introduction 
of temperature boundary conditions. In addition we added an external forc- 
ing term and removed the Kepplerian shear and the accompanying it shearing 
coordinates. 

The resulting code is ideal for our purposes for the following reasons: 
1. The anelastic approximation means no shock or sound waves can exist 
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and hence the time steps can be much larger than for a fully compress- 
ible code, and yet the anelastic approximation allows phenomena such 
as convection to occur unlike in completely incompressible codes. This 
allows us to get much closer to the actual time scales of interest (usu- 
ally of order hours or days) than a fully compressible code would. The 
Boussinesq approximation has those same properties, however it is not 
suitable for simulating stratified convection where the buoyancy does 
not depend entirely on the temperature. 

2. Spectral codes are more efficient at simulating turbulent flows than flnite 
difference codes because their spatial accuracy is exponential rather than 
a power law, so they are capable of reliably reproducing the turbulent 
flow even at modest resolution. 

3. The grid happens to have higher density near the top and bottom of 
the box where it is needed in order to resolve the boundary layers that 
develop there due to the boundary conditions we impose. 

In this paper we present the details of the modifled code as well as a 
perturbative estimate of the efficiency of turbulent viscosity. In section [2] we 
derive the anelastic approximation to the Navier-Stokes equations modified 
to include a time dependent shear. In section [3] we present the numeri- 
cal algorithm. In section H] we show the tests we ran to confirm that our 
implementation actually evolves the desired equations, and verify that the 
numerical scheme is indeed second-order accurate in time. Finally, in sec- 
tion O we use a perturbative calculation to get an estimate of the effective 
turbulent viscosity in our simulated box. 

2. The evolution equations 

We define upfront the following quantities describing the fiuid and the 
fiow: 



p 


- pressure 


T 


- temperature 


p 


- density 


V 


- velocity 


e 


- potential temperature 


Cp 


- constant pressure specific heat 


R 


- ideal gas constant 


9 


- acceleration of gravity 




- heat diffusion coefficient 


f 


- some external forcing (e.g. tidal force) 



As was said before we do not actually evolve the fully compress i ble fiu id 



equations but rather their anelastic approximation (c.f. iBannonI (119961 )) 
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Under the anelastic assumption each flow variable is split into a sum of a 
time-indcpcndcnt background component (denoted by an over bar) and a 
time varying perturbation (denoted by tilde over its symbol) superimposed 
on top of that. 

The background variables satisfy the fully compressible equations with all 
X, y and t derivatives set to zero. The background velocity is assumed to be 
zero as well. The boundary conditions that complete these equations are that 
we require the background temperature on the top and bottom walls to have 
some fixed values Ti and Th respectively, and the pressure at the top wall to 
have some value, ptop- Alternatively, we could fix the temperature gradients 
at the boundary, which would set the flux that the simulated convective 
layer must transport. Neither of those two options is exactly what the actual 
physical problem requires. We could base our boundary conditions or reliable 
physics if we were planning to simulate the entire convective zone. However, 
this is impossible to do at the resolution required to study the turbulent 
dissipation, so we will ultimately be interested only in the flow that develops 
in the interior of the box and plan to exclude the regions near the top and 
bottom boundary from any analysis we do. For this reason the choice of the 
exact thermal boundary conditions will hopefully have little impact on our 
results. 

For convenience, we define the following quantities: 



Where the z coordinate has a value —Lz/2 at the bottom of the box and 
Lz/2 at the top of the box. 

The background variables that emerge as the solution to fully compress- 
ible equations are: 




(3) 
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P 



pRT 



(7) 



R 




(8) 
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With those definitions the anelastic equations governing the evolution of the 
perturbation variables become: 



V-pv = 





dv 
'dt 


V X c<j — 


ae 

di 


de 

^ dz 


h = 


p V 


p 
p 


P T 


e 


p 


e 


pgHp 



(9) 

Vh + + f (10) 
~~ ^ :V ■ (kVt) (11) 



(12) 
(13) 
(14) 



Where we have introduced an enthalpy h (defined by equation [12]) , a density 
scale height Hp = — (^^^) and the vorticity = V x v. Equation [H] is 
not the correct linearized equation for the potential temperature. The de- 
nominator of the pressure term should have be en 'yp instead of pgHp. This 



replacement was introduced by iBannoru (119961 ) . He showed that this sub 



stitution is required to ensure that the anelastic equations conserve energy. 
To define the time evolution completely we need to add boundary conditions 
to the above equations. The boundary conditions on the four vertical walls 
of the domain are set by the fact that we use Fourier series expansion for 
the horizontal spatial dependence of all quantities, and hence all quantities 
are naturally periodic in those directions. In addition to that, we want the 
temperature at the top and bottom boundary to be whatever we specify, 
and hence its perturbation should be zero. Also we require that the top and 
bottom walls are impermeable, but with no friction. That means that we set 
Vz to zero at the top and bottom boundary, and do not require anything for 
Vx and Vy. 
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The anelastic equations above obey a set of energy conservation equations 
for the following definitions of the kinetic and thermal energies: 

Ek = [ p^dV (15) 



Ej 



1^ C.-pfidV 



(16) 



Using the anelastic evolution equations we can express the rates of change of 
these two energies to be: 

dEK 



dt 

dEx 
"dT 

Where sinks/sources are defined as: 
£i = g 

£2 = 

3. Numerical Time Evolution 



£1 



-£i + £2 




dT 

dz 



dxdy 



(17) 
(18) 

(19) 
(20) 
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3.1. Spectral Method 

The spectral representation for a flow variable (q) used in the code is 
given by: 



Nx 
2 



Ny 
2 



q{x,y,z,t) 



J2 E 5^gfc^n^(t)e^'"'^/^^e^^-'^/^^T„ 



(21) 



Where the vertical basis functions: Tm(z) = cos ( m cos ^ 



, are Chebyshev 

polynomials. The spectral method we use does almost all calculations in the 
wavenumber/ Chebyshev space, except for taking products of variables, which 
are done by first transforming back to physical (x, y, z) space on a grid of 
collocation points, taking the product there and trans forming back. For a 
more complete discussion of the spectral method used see lBarranco &: Marcus 
f l2006h (section 3 and 3.2). 
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3.2. Time integration 
3.2.1. Advection Step 

The advection step is fully explicit. It uses variables from this and the 
previous time step to achieve second order accuracy. It first calculates the 
quantities: 



9 

Tl = V X LU + -gz 
9 

^ = -v,^-^-V9 

dz 

Then velocity and temperature are updated using: 
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N+i 



2 

2 ^ 



N 



(22) 
(23) 

(24) 
(25) 



3.2.2. Hyperviscosity Step 

The purpose of this step is to suppress the highest modes both in the 
horizontal and vertical directions to avoid buildup of power there due to the 
truncation of the spectrum at some finite number of spectral coefficients. For 
finite difference codes a step of this sort is unnecessary because there is some 
finite "grid viscosity" associated with the numerical scheme. In spectral 
codes there is no equivalent effect that prevents the build up of power at 
the highest simulated wavenumber r aodes. We implemeri t the h ypervisocisty 
step exactly in the way described in lBarranco &: Marcus! (120061 ) . suppressing 
each spectral coefficient by a factor as follows: 



N+- 

V 3 



qN+\ 



N+- 

V ^3 exp 



Q^+\ exp 



-At ( v'^yPk'^ + i/^'^Pm^P 



(26) 
(27) 



Where and v^^^ are some hyperviscosity coefficients, p is an integer 
between 1 and 6. k\ = + k^ is the horizontal wavenumber and n is the 
order of the Chebyshev polynomial that this particular amplitude applies to. 

3.2.3. Pressure Step 

The pressure step basically ensures that the velocity at the end of the time 
step satisfies the anelastic constraint. In the numerical scheme we abandon 
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the enthalpy at a specific time as a variable, and instead use its average 
between two consecutive time steps: 



(28) 



We then achieve second order accurate time evolution by updating the ve- 
locity according to: 



At 



N 



N+l 



(29) 



So imposing the anelastic constraint we get: 



= v-v^+^+.r^^ 

dz 



dz 



dz dz 
(30) 



Regrouping and imposing f^L lz = we get the differential equation for 

^ 2 

updating 11 with its boundary conditions: 



^2 rflogp d 
V + 



dz dz 

QllN+l 



n 



N+l 



At 



N+^dlogp 
dz 



dz 



1 Af+I 
Vz 

At 



±L^/2 



(31) 
(32) 



For details of the implementation of this equation in our code see appendix 



3.2.4- Heat Diffusion Step 

The heat diffusion step updates the potential temperature according to: 



e 



'N+l 



At 9Rn 



d\\in d 

dz dz 



rpN _j_ rp 



N+l 



(33) 
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When we express T in terms of v,/i and 6, this gives a second order differential 
equation for 6^~^^. In terms of the following definitions: 





- - 


gnf 




dlnn 


2g 


dz 




g f dlnO 




dz 


2Cpp 





dz dz 



dlnK, 



+ 



K 



K 
K 



h- 



dz 



K,At 



gnT 



(34) 
(35) 
(36) 
(37) 
(38) 



The equation we solve during this step is: 



(<p^+i+<p^) 



d_ 
dz 



C2 1 9^ + c.e 



(39) 



This is a second order equation for 6^^^ so we need two boundary condi- 
tions to make the solution unique. These come from the requirement that 

the temperature perturbation on the boundary vanishes: T . = 0. This 
requires: 



e 



N+l 



C, ( - - 



(40) 



So we see that the differential equation [39] uses only the value of 11^+^, but 
the boundary conditions need the value of the enthalpy at the updated time. 
There are two options of how to obtain this value. One is to keep track of 
the enthalpy from the very beginning and after each pressure step to update 
it as: ^ ^ 

^N+i ^ 2n^+i _ Z,^ (41) 

However, if the initial value of h is not perfectly set this prescription will lead 
to oscillations in the value of the enthalpy at the top and bottom boundaries. 
To illustrate this assume that initially we set to a value that is slightly 
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higher than what it should be. Since 11^ is calculated without reference 
to the initial conditions it will have the correct value. This will lead to 
being slightly smaller and so on. These oscillations are then translated to the 
interior of the box through their effect on the potential temperature variable. 
Further, we found that these perturbations grow with time for all the cases 
we ran. 

As a result instead of introducing the additional variable h with the only 
purpose of getting temperature boundary conditions we use an approximation 
to its value through 11: 

^ ^ (3n^+^ - n^) (42) 

These boundary conditions are then imposed by ignoring the differential 
equation at the top and bottom wall (where it does not make sense any way) 
and computing the values of 9 at those boundaries to satisfy the boundary 
conditions. For the implementation details of this step see appendix [B] 



4. Tests 



Since all the Fourier trans forms and differentia l opera tors are computed 
in exactly the same way as in lBarranco fc MarcusI ( l2006l ). see section 4.1 of 
that paper for a discussion of the performance of the code. 

Below we present various tests we ran to confirm that the code is indeed 
evolving the equations described above, and that the numerical scheme em- 
ployed is indeed second order accurate in time. We repeat essentially all the 
tests that Barranco and Marcus ran to confirm their numerical scheme. 



4.1. The background state 

Most of the tests we performed have the same background state, which is 
also the state we used for the final convective box simulation. The only time 
we modified the background was in order to simulate a convectively stable 
box and observe g-modes. Here we describe that background state. We use 
a self consistent non-dimensional set of units. Hence, the results apply for 
any choice of units. 

The dimensions of the box are 4x4x4, with a resolution of 128 collocation 
points per direction. 

Next we need to choose the vertical profile of the heat diffusion coeffi- 
cient. On one hand, because we require the vertical velocity to vanish at the 
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Figure 1: a) The heat diffusion coefficient (k), b) the background potential temperature 
{6), c) the background density (p) and d) the background pressure (p) profiles. 



top and bottom boundaries, we need to set up a background that is stable 
to convection near those boundaries. On the other hand, we would like to 
simulate a top-driven convection, which is typical for stars with surface con- 
vection zones. To achieve that we need the most unstable stratification to 
occur near the top of the box. Further, since the heat diffusion step uses 
the heat diffusion coefficient, and its first and second derivatives, we need 
our expression for the heat diffusion coefficient to have a continuous second 
derivative. 

To achieve all these requirements we construct the particular k,{z) profile 
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used in this paper from 6 separate pieces as follows: 



Ki + K2 Sm{k{z - Zq)) 



,Z < Zq 
,Zq< Z < Zi 



,Zi < Z < Z2 



(43) 



v 




Z2)^,Z2 < z < z^ 



,Z'i<Z < 2:4 



,^4 < Z 



Where the parameters Ki, i = 1 ... 12, Zj , j = 1 ... 4 and k are chosen to make 
K, k' and k" continuous and to select the shape of the curve. The shape of 
the curve used in this article was determined from the following constraints: 

• The values that n{z) takes at the boundaries: = 20 and = 21. 

• The depth above which k{z) remains at its maximum value of Ki2'- 
Zi = 1.95. 

• The depth at which k,{z) has a minimum and the value at that mini- 
mum: k{z2 = 1.4) = 19.8 

• The depth below which k,{z) is held constant at its bottom value of kq: 
zo = -1.8 

• The locations of the two inflection points: zi = 0.2 and z^ = 1.85. 

A plot of the resulting depth dependence of n{z) is presented in figure [1^. 

Next we need to choose values for the background temperature at the top 
(Tiow) and bottom (Thigh) of the box. Those are dictated by the requirement 
that the flow speeds should never approach the local speed of sound; other- 
wise, the anelastic approximation is no longer an acceptable approximation 
and fully compressible equations should be used. 

Finally we need to choose values for the pressure at the top of the box 
(Ptop), the external gravity (g), the specific heat at constant pressure (Cp), 
and the value of the ideal gas constant (R) . These values, together with Tiow 
determine the pressure and density scale heights. Since we are interested in 
studying turbulence in a stratified medium we need to choose these values 
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such that our box encompasses at least several pressure and density scale 
heights. The particular values we chose were: 



Ptop — 


1.0 X 10^ 


(44) 


9 = 


2.74 


(45) 


Cp = 


0.21 


(46) 


R = 


8.317 X 10^2 


(47) 


Tlow 


10.0 


(48) 


Thigh 


62.37 


(49) 



We deliberately do not include any units in the above quantities, since the 
simulated flow is independent of the choice of units. 

The resulting background potential temperature, pressure and density are 
presented in figures [DDcd. One can see that the steepest negative slope of 9 
occurs for heights between 1 and 1.5 units. Also we see we have significant 
stratification. From figures [D^d we see that the convective box encompasses 
2.8 density scale heights and 4.6 pressure scale heights. 

J1..2. Energy 

We first verify that the energy like conserved quantities defined in equa- 
tion [15] and equation [16] evolve according to eqations [T7] - [20] The initial 
conditions we used for this test were that all components of the velocity were 
set to 0, and the initial potential temperature contained random fiuctuations 
in the lower 10% of the spectral modes. We chose the time step to be much 
smaller than the smallest absolute value of the buoyancy period, which in 
this case is an imaginary quantity for most of the box. This way the time 
step is short both compared to the growth rate of the instability near the 
middle of the box and the g-mode period near the boundaries. 

We perform the test by running the code with no hyperviscosity. We 
output the kinetic and thermal energies as well as the rates £1 and £2 at 
every time step. We then compute the time integrals of £1 and £2 using 
Simpsons's method. This is sufficient since the numerical evolution is only 
second order accurate in time. In figure [2] we show that equations [T7] - [20] 
are indeed satisfied to one part in a million. 
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Figure 2: Energy conservation. Left: negative thermal energy - solid line, the absolute 
difference between the total energy calculated directly and as the integral of £2 - dashed 
line; Right : Kinetic energy (solid line), the absolute difference between the Kinetic energy 
calculated directly and as the integral of £1 - dashed line 



4.3. Normal modes 

Another test we ran was checking that the normal modes of the hnearized 
equations of motion evolve as expected. We look for normal modes of the 
form: 

y, z, t) = q{K, ky, ^)e— (50) 

For convenience define 



d 

T = 



<^B = 9 



dlogO _ 9_f 9_ _ Thigh - Tio^ 

dz ~ f\Cr, 



(51) 

(52) 
(53) 



In terms of equation [5T] and equation [52], and after dropping all nonlinear 
terms, the anelastic equations (P- [T^ simplify to include only two variables 
V. and f: 



— iuj { I 



-lUJ 



DD, 



f+^{D^-D)v, 
Cpp 



K 



Cpp 



-kl + iD + D^ 



(54) 
f55) 



Where the operators D and Da are defined in appendix A equations [73] and 
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z z 

Figure 3: The shape of the normal modes we initilized the box with. On the left is the z 
component of velocity and on the right is the potential temperature. 



In terms of Vz and f the remaining flow variables can be expressed as: 

iujDaVz (56) 
ikji (57) 
ikyh (58) 

We calculated numerical solutions to the eigenmode equations (see figure [3]) 
and supplied them as initial conditions with very low amplitude to the code. 
We expect that for later times the evolution is done simply by multiplying 
the initial amplitudes by e~*'^*. We ran the code with the background already 
discussed above and a time step that was no larger than for the mode in 
question. The comparison with the simulated evolution of these eigenmodes 
is presented in figure |H As we can see for large enough time steps the error 
scales as the square of the time step, which confirms that the code is indeed 
second order accurate in time. For very small time steps the error deviates 
from that scaling due to numerical round off. The minimal fractional error is 
much larger than the numerical precision, because we have chosen the heat 
diffusion to be as small as possible, which causes the matrix we invert during 
the heat diffusion step to have values along the diagonal that are many orders 
of magnitude larger than the off diagonal values, which causes the numerical 
roundoff to be amplified many times. This also explains why the error in 
the potential temperature is largest (the other varibles suffer from this only 
indirectly). 



klh = 
iojVx = 
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Figure 4: The maximum error in normal mode evolution, after a fixed evolution time. 
The error is scaled by the maximum value of the quantity at time t = 0. The solid line 
corresponds to quadratic scaling of the error with the time step. 



4-4- 2D vortices 

Another test we ran was initializing the box with a pair of circular columns 
of vorticity running from the top to the bottom boundary. In physical space 
for each column the vorticity was constant and in the z direction. The two 
columns had opposite signs of vorticity so as to ensure that the total vortic- 
ity in the box is zero, as required by the periodic boundary conditions. One 
expects that the two vortices will move parallel at a constant rate determined 
by the distance between them and the magnitude of their vorticity. Because 
of the periodic boundary conditions, having two vortices in our box is equiv- 
alent to having an infinite number of vortices, copies of the two in the box. 
There is no analytical expression of the infinite series for the velocity with 
which each vortex should move, but if the two vortices are far away from the 
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20 40 60 80 100 

time 

Figure 5: The difference between expected vortex position and the average position of 
vorticity. The period of the osciUations corresponds to the rotational period of the vortices 
and the drop in the end is due to the fact that the vortices are exiting the box on the right 
and hence a bit of them is appearing in the left, causing the average position of vorticity 
to move toward zero. 



walls as compared to the distance between them we expect that their motion 
is at least approximately that of the situation of only 2 vortices. The time 
step we chose in this case was O.Oltcross (the expected time it would take the 
vortices to cross half the box). We ran this test and verified that the rate at 
which the vortices moved corresponded to the approximate analytical rate, 
see figure |3 

4.4- 1- Convectively stable box 

The last test we ran was imposing a convectively stable stratification in 
the box, and initializing with random temperature perturbations. In this 
case one expects to see g-mode oscillations with a frequency given by the 
buoyancy frequency cu^ = gd In 9/ dz. So for this test we needed the buoyancy 
frequency to be approximately constant throughout the box, and much larger 
than the time step, but small enough to allow us to simulate many buoyancy 
periods. So for this test we set the heat diffusion coefficient to be constant 
throughout the box. Also we set the top and bottom temperatures to be 
the same. This way the buoyancy frequency was independent of height. The 
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Figure 6: Kinetic energy as a function of time for a convectively stable box. The distance 
between every two consecutive vertical lines is the buoyancy period. The time is in units 
where the buoyancy period is 1. 

particular values for the parameters of this run were: 



Ptop = 


1.0 X 10^ 


(59) 


9 = 


2.0 


(60) 


Cp = 


0.2 


(61) 


R = 


8.317e-2 


(62) 


'-p 

low 


10 


(63) 


Thigh 


10 


(64) 


dt = 


0.002 


(65) 



1 



This means that the buoyancy frequency is u;^ = 2.0, A plot of the kinetic 
energy for this run is presented in figure O It can be seen that the kinetic 
energy goes through two cycles every period. This is because the velocity 
has to go through one cycle, and the energy has a maximum every time the 
velocity has a maximum or a minimum. The decay in the curve is caused 
mostly by the heat diffusion smoothing out the perturbations over time. 
Since different modes decay at different rates and they are coupled through 
the nonlinear terms we have not shown an expected decay curve. 
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5. Estimating the dissipation in an unstably stratified convective 
box 



5.1. Steady state flow 

Having confirmed that the code is solving the correct equations, we ran 
a box with the background state described in section 14.11 for long enough to 
reach a steady-state flow. The criteria for having reached a steady-state flow 
were that the kinetic and thermal energies should stop drifting systematically 
up or down (see figure [TT]). The oscillations we see in the kinetic energy have 
a period close to the convective turnover time of the box. We also want 
the spatial spectra of the velocity and potential temperature to be steady to 
within a few percent. The vertical spectra can not be directly obtained from 
the output of the simulation since the collocation points of computational 
grid are not evenly spaced in the vertical direction. So we first had to re- 
sample to an evenly spaced grid and apply some window function in the 
vertical direction before performing the discreet Fourier transform. Since in 
the vertical direction we simulate the Chebyshev series of the quantities the 
most natural way to re-sample to an even grid was to evaluate this series for 
each of the new (evenly spaced) collocation points. Fourier power spectra 
of the 3 velocity components and the potential temperature are presented in 
figure [TJ along with the scaling that Kolmogorov statistics predict (spectral 
power oc k~^^^). The sharp cutoff at high wavenumbers is related to the 
resolution of the box (using a box of half the resolution produces a cutoff at 
half the wavenumber we see in figure [7]). 

Horizontal cross sections of Vz and ^ + ^ at 3 different heights of a typical 
steady state frame of the flow are presented in flgure [3 The three heights 
we chose were z = 1.98, z = and z = —1.98 for a box in which the vertical 
coordinate runs from -2 to 2. Horizontal section of the z component of the 
vorticity Uz as well as typical vertical sections of the above quantities are 
presented in flgure M 

At the top of the box (flgures [Hfe',d), near the top boundary, the flow 
exhibits a cellular pattern of narrow downflow lanes. Traces of the cells are 
visible only in the top 5% of the box, after that the flow organizes itself 
into a pair of perpendicular downflow lanes that horizontally span the entire 
box. With time, those lanes get distorted, break up and reform, but are well 
deflned for at least half the time, for most of the upper half of the box. They 
ar e generally parallel to the g rid axis. Similar patterns were flrst observed 




Their tests show that the lanes tend to align 



20 




Figure 7: The (x: top left, y: top right, z: bottom left, time: bottom right) spectra of the 
3 velocity components and the potential temperature. The thick line in the spatial spectra 
plots corresponds to Kolmogorov scaling {Ek oc k~^^^). The thick line in the time spectra 
plot corresponds to the scaling we find for the effective viscosity in the next section. 
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Figure 8: Horizontal cross sections of the vertical velocity component (left) and the total 
potential temperature (right) at three different heights of the box: z — 1.98 (a,d), z = 
(b,e) and z = -1.98 (c,f) 
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Figure 9: Horizontal cross sections of the vertical vorticity component (left panels) at 
heights: z = 1.98 (a), z = (b) and z = —1.98 (c). Vertical cross sections (right panels) 
of the vertical velocity (d), the total potential temperature (e) and the vertical vorticity 
(f)- 
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Figure 10: The average logarithmic gradient of the potential temperature. 



themselves with the periodic directions of the models, and not with the grid 
axis per se (rotating the periodic directions at 45° relative to the grid axis 
caused the lanes to rotate as well). So they conclude this to be due to the 
small horizontal to vertical aspect ratio of the simulations. 

Around the middle of the box the pair of perpendicular lanes are still 
sometimes visible, but they are a lot less persistent. Rather at those depths 
the flow consists of one or two downflows, which take up less than a quarter 
of the cross sectional area among a gentler upflow (figures [8]3,e). For the 
lower half of the box, the asymmetry between up and downflows gradually 
decreases as we get further down. Near the bottom boundary the large scale 
of the flow disappears, until only small scale structure is left near the bottom 
boundary (figures [Ht,f). 

Since we ran our simulations with the smallest possible value of the heat 
diffusion coefficient we expect to have very efficient convection, hence we 
would expect the z gradient of the potential temperature to be very close to 
zero except near the impenetrable top and bottom boundaries. The averaged 
over time and horizontal directions logarithmic gradient oi 6 + 6 can be seen 
in figure [TOl As we can see the scale height of 6 is indeed more than two 
orders of magnitude larger than the box, as long as we are not very close 
to the boundaries. The top of the box has a local minimum of the entropy 
gradient and a significant convectively stable region, so as to keep the large 
magnitude vertical fiow from getting too close to the boundary and causing 
numerical problems. We also see that another local minimum of the entropy 
gradient has developed near the bottom of the box. This is caused by the 
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Figure 11: Kinetic (left) and thermal (right) energy content of the convective box used as 
one of the criteria for having reached a steady state. We decided steady state was reached 
at time time of 400. 



lower impenetrable wall. Not having a significantly thick convectively stable 
region near the bottom of the box is acceptable, because the larger density 
of the fluid means that the flow velocities are much smaller at the bottom 
than at the top of the box. 

5.2. Lowest order perturbative expansion 

Ignoring the part of th e run before steady stat e was reached we adapted 



the lGoodman fc OhI fll997h method as described in lPenev et aP fl2007l . l2008bl ) 
to find a lowest order perturbative expansion estimate of the energy transfer 
rate between a small amplitude external forcing and the turbulent flow in 
our box, and respectively from that we can derive the components of an 
effective viscosity tensor that reproduces the energy dissipation rate due to 
the turbulence. We assume forcing in the form of an external velocity field: 

V = A{t) ■ X (66) 

Goodman fc OhI ( 1997 ) define two dimensionless parameters: the tidal strain 
|A|, and (fir^)"^, where Q is the frequency of the perturbation and 
Tc = Lc/Vc- The characteristic convective length scale is and Vc is the 
characteristic convective velocity. In the case of hierarchical eddy structured 
convection Tc is the largest eddy turnover time. 



We then take only the lowest nonzero term in the expansion of the secular 
rate of change of the kinetic energy £ with respect to those two parameters. 
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(jPenev et al.l . l2008bl ) adapt the (iGoodman fc OhI . 119971 ) formalism for dis- 
cretely sampled data on a finite spatial and temporal domain the resulting 
expression is: 



2T r r ^ 



-R 



X,IJ,,u,u' 
4 1 



. ^-"^ CJn I. L J J 



X,IJ,,u,u' PT^O 

Where the simulations are assumed to span a time in the range of (0, T); N^, 
Ny, Nz and Nf are number of grid points in the x, y, z and time directions 
respectively; Af = N^NyN^Nt, A, /i, z/, p are indices for grid points in discrete 
Fourier space: k^^^,;/ = {2n\/ Lx.lnp/ Ly,2T[v / Lz), Wp = 27ip/T; v is the 
time and space Discrete Fourier Transform (DFT) of the convective velocity 
in the absence of the external perturbation (i.e. the simulated velocity field); 
p is the DFT of the background density along the z direction; Px^^^u = 
I — \i.\^i,^,M\^ii,^,Jk^_^^ j, is the discrete version of the Pk operator defined by 



( iGoodman fc OhI . 119971 ) that maintains incompressibility. 

The 2, 2 subscripts of £'2,2 denote that equation [67] contains up to second 
order terms in the two dimensionless parameters characterizing the tide and 
the convection respectively, p{kz) is the Fourier transform of the density 
averaged over x, y, t. The normalization is such that p(0) is the average 
density over all space and time. 

To perform the calculation we use a discreet Fourier transform to estimate 
the spectra of the velocities and density needed in equation |67l but before 
that we again re-sample to a grid that has its collocation points evenly spaced 
in the vertical direction. We have to be careful when using discreet Fourier 
transforms to estimate spectra. In particular we need to pay special atten- 
tion to the z and time directions, since they are not periodic and hence the 
discreet Fourier transforms suffer from windowing effects. To confirm that 
our results are not significantly affected, we perform the same calculations 
using no windowing in those two directions, as well as Welch (parabolic) and 
Bartlett (triangular) windows. Also it is possible that the impenetrable top 
and bottom walls might influence our result. So we also tried a Welch and 
Bartlett window functions that exclude the top 12.5% and the bottom 5%. 
The resulting effective viscosity scalings from all these tries were not signifi- 
cantly different. 
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Figure 12: The four in dependent components of the viscosity tensor as estimated by the 
Goodman &: OhI (1997) method. The vertical axis is the effective viscosity in units of 



< vt 



>^/^ Hp {Hp - pressure scale height) and the horizontal axis is the perturbation 



frequency is units of inverse convective turnover time (r) 



We then note that the energy transfer rate of equation [67| is due to the 
external forcing. If the flow is steady state this transfer rate must be balanced 
by the energy dissipated by the turbulence. We can then extract the different 
components of the effective viscosity by setting all terms of A to zero except 
for one, and matching to the energy dissipation rate that an actual molecular 
viscosity {u) would cause: 

S^isc = \ {pv) Trace [A(n) ■ A*(fi)] , (68) 

where averaging is done over the volume of the box and over time. Note that 
this is only the dissipation of the flow caused by the external shear. In the 
absence of external shear (i.e. A — 0) there is no dissipated energy. The 
dissipation of the convective flow is not of interest because it is assumed to 
be balanced by the thermal driving of the convection. 

Since physically there should be no difference between the two horizontal 
directions we expect that v^x ~ ^yy, ^xy ~ i^yx and v^z ~ i^yz ~ ^zx ~ i^zy 
We verified that these approximations indeed hold. So we are left with only 
four independent components of the viscosity tensor. Those are presented 
in figure [121 As we can see, all the components have similar scaling with 
frequency, and the diagonal components are a bit over an order of magnitude 
larger than the off-diagonal components. The approximate frequency scaling 
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for the low frequency dependence of the viscosity components (for Q < 100) 



is: 



u oc n'-'^'-' (69) 



6. Conclusion 



The above result points to the possibility that viscosity in turbulent con- 
vection zones loses efficiency significantly slower than what Kolmogorov scal- 
ing predicts at least on large timescales. This seems to be due to the fact 
that the turbulent eddies with turnover times similar to these large timescales 
are not in the inertial subrange and hence, the velocity power spectrum is 
much shallower than the Kolmogorov 5/3 law. This result is not conclusive, 
since the possibility remains that dropping higher order terms in the above 
expan sion is not a good ap proximation. 



In iPenev et al.l ( l2008af ) we introduce horizontal depth dependent forcing 
into the fiow equations and obtain the dissipation properties of the turbulent 
convective fiow directly. We were then able to compare the average rate of 
work done on the fiow by the external forcing to the expected rate of energy 
transport and dissipation by an assumed effective viscosity. We found that 
with sufficiently long time average these two quantities have the same depth 
dependence, thus verifying the validity of the assumption that an effective 
viscosity coefficient is sufficient to parametrize the average dissipation and 
momentum transport properties of the turbulent convective fiow. Further, 
by repeating the above procedure a number of times we were able to derive 
the scaling of this effective viscosity coefficient with period and confirm that 
it is linear. 

In this paper we have presented a code that is well suited for the purpose 
of studying turbulent dissipation in convective zones. First it is a spectral 
code, which means that the spatial accuracy is exponential, and hence the 
code efficient for simulating turbulent fiows. Further by using the anelastic 
approximation we do not need to deal with sound waves and shocks. This 
allows us to take much larger time steps (by more than an order of magnitude) 
than with a fully compressible code, and hence we can have runs that cover 
a much longer physical time than fully compressible simulations. This is 
important since the external forcing regimes we are interested in are tides 
with orbital periods of several days, which is not achievable in reasonable 
time with fully compressible codes. 
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The price we pay of course is that the flow we simulate is not a good 
approximation to the flow in any star's surface convection zone. Firstly, we 
cannot accommodate the region where the driving of the convection occurs, 
because this region is characterized by strongly supersonic flows and shocks, 
also the mean free path of light in that region is not small and hence the 
radiative effects can no longer be captured by a heat diffusion coefficient. 
Secondly, our code uses the ideal gas equations of state which is a poor 
approximation to the upper layers of stars. 

However, the scaling of the effective viscosity with frequency we obtained 
(equation [69]) for our box agrees with the scalings obtained with the same 
perturbative approach applied to realistic simulations of the upper layers of 



the co nvection zones of the sun and two smaller stars (jPenev et al.l . 120071 . 



2008bl ). Which gives us confidence that our results are applicable to the 



systems we are interested in studying. 

We would like to acknowledge the contribution of Philip Marcus to the 
development of the original code (IBarranco fc Marcusl . l2006l ). 



A. Implementing the Pressure Step 



As discussed in section 13.2.31 in order to advance the enthalpy by one 
time step we need to solve the following differential equation with boundary 
conditions: 



dlogp d 

dz dz 



N+l 



dz 



n 



-L ( v-v 

At 
1 

At 



+ V 



N+l dlogp 
dz 



(70) 



±L,/2 



The solution is obtained in two steps. First we impose the boundary con- 
ditions by ignoring the differential equation for the two highest Chebyshev 
modes and replacing it with the equations for the boundary conditions. Then 
we use a Green's step (see sec. lA.ip to fix the equation for those two high- 
est modes and instead impose the boundary conditions at the expense of 
satisfying the equation at the physical top and bottom boundary. 

We need to solve this equation as efficiently as possible. Clearly we can 
make this a trivial matrix multiplication operation if we were to store the 
inverse of the left hand side operator for each horizontal mode at the start, 
and at each time step we multiply every vertical slice of the transformed right 



29 



hand side by the corresponding inverse to get the value of n^+^. However, 
this would require a A^^. x A^^^ matrices of N"^ elements to be stored (where 
Ny and A^^ are the resolutions in the x, y and z directions respectively), 
which for large resolutions is likely to exceed the amount of memory available 
on each node of the cluster where the code is to run. 

To avoid this we need to find the most efficient way which only stores 
things common to all horizontal modes. We see that the pressure equation is 
almost identical for all horizontal modes, except for the horizontal component 
of the operator, which in Fourier space means simply multiplying by 
k']_ = kl. + ky, where kx and ky are the corresponding x and y wavenumbers 
for each mode. So what we can reasonably do is pre-compute the common 
part of the left hand side operator and then for each horizontal mode add k^ 
along the diagonal. We then overwrite the last row of the resulting matrix 
(Mi,j) with 

Mn.,, = (71) 
M^,,_i,^ = (-l)V (72) 

in order to impose the boundary conditions. We then decompose M into an 
upper and lower triangular parts and solve by backward substitution. The 
right hand side vector also needs to have its highest two entries overwritten 
with the boundary conditions at the top and bottom of the box respectively. 

A.l. Green's Step 

In this step we fix the anclastic constraint even for the two highest Cheby- 
shev modes and instead break the two highest modes of the momentum equa- 
tion to satisfy the boundary conditions which we break in the process. To 
make expressions shorter define the following operators: 



D = 


d 
dz 


(73) 




d dhip 


(74) 


Da = 


dz dz 


V = 


ikx^ + ikyY + Dz 


(75) 


Va = 


ikx^ + ikyY + Daz 


(76) 


Aa = 




(77) 



We would like to modify the pressure step in a way that will include two new 
degrees of freedom which we can then use to fix the anelastic constraint for 
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the two highest Chebyshev modes. The particular modification useful in this 
case is: 



^N+i ^ ^N+l _ AiVn^+^ + z (rf +^Tm-i + t^^'Tm) (78) 

Where T/v denotes a Chebyshev polynomial of order A^, ri and T2 are arbitrary 
coefficients to be chosen later and M is the order of the highest Chebyshev 
coefficient we are simulating. 

Requiring the anelastic constraint and velocity boundary conditions gives: 

A^n^+^ = -^WA-y''+'^+r,^^'DATM-i + r^+'DATM (79) 



To proceed we break up n-^+^ into 3 pieces: 



±L,/2 



(80) 



niv+i ^ n^+i ^ 1 (^^^+irf +1 + r^+^r^+^) (8i) 

This allows us to split the above equation into 3 separate equations with cor- 
responding boundary conditions, the first of which is the already calculated 
pressure step: 

AaK"-' = ^V^-v^+i (82) 



tN+11 _ _L„^+I 

^ l±W2 - A^^- 



(83) 

±L./2 

A^rf+i = DaTm-1 (84) 

D^l^^\±L,/2 ^ ^M-l|±L^/2 (85) 

A^r^+1 = L>^rM (86) 

^^2'^^\±L./2 = ^m|±l^/2 (87) 

We solve the two new equations the same way we solve the first. This makes 
as before the anelastic constraint satisfied for all but the two highest Cheby- 
shev modes, but this time we have 2 arbitrary constants Ti and T2 which we 
can then set to values that will make the anelastic constraint hold for those 
two modes as well. 
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B. Implementing the Heat Diffusion Step 

The same considerations as the pressure step apply to this step. We 
again decide against making a table of pre-inverted matrices for each hori- 
zontal mode, and instead we only pre-compute the part of the matrix that is 
common for all horizontal modes. To make the description of the numerical 
procedure followed we define the following matrices: 
CP:transforms a vector from Chebyshev to physical space 
PC:transforms a vector from physical to Chebyshev space 
D : differentiates a vector in Chebyshev space 
From those we construct a derivative operator and the common part of the 
left hand side operator, both in physical space: 

Dp = CP-D-PC (88) 
Mp = Dp- Dp + F ■ Dp + G (89) 

What we pre-compute and store is a matrix M, which we obtain by 
overwriting the first and last row of Mp with: 

Mpi,i = 1 Mpi,, = 0, J = 2..N, (90) 
MPN.,, = 0, J = l..iV, - 1 Mpjv,,jv. = 1 (91) 

to allow for imposing the boundary conditions, then sandwiching the result- 
ing matrix between PC and CP: 

M = PC -Mp- CP (92) 

Then in order to solve equations [31] and SD] for each horizontal mode we 
construct a new matrix 'WL'{k±) from M as follows: 

M:^^.(A;x) = M,,, - kl6,,, + kl {PC,,oCPo,j + P^.^v.-iCP^.-i,,) (93) 

The first new term adds the horizontal part of the operator for the given 
mode. However, this breaks the requirements for the boundary conditions, 
so we repair them with the other two terms. This way the boundary con- 
ditions are directly imposed by breaking equation [31] at the physical top 
and bottom of the box, instead of ignoring it for the two highest Chebyshev 
modes. That means we do not need to perform extra greens steps like we 
did for the pressure equation. Then we obtain the solution to equation [391 by 
decomposing M'(fcx) into an upper and lower triangular matrices and using 
back-substitution. 
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